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1.  INTRODUCTION 


In  composite  laminates  made  of  high  stiffness  and  strength  fiber-reinforced  plies,  the  defor¬ 
mation  effects  due  to  transverse  shearing  and  transverse  normal  stretching  can  often  be  signifi¬ 
cant,  especially  in  thick  laminates  and  those  subject  to  short  wavelength  loading.  These 
effects  are  more  pronounced  in  composites  than  in  homogeneous  materials  due  to  their  inher¬ 
ently  high  material  compliancy  in  the  transverse  directions  relative  to  the  axial  fiber  direction. 
Furthermore,  composite  laminates  exhibit  much  lower  strength  in  the  transverse  directions  and 
ar  the  ply  interfaces,  thus  being  particularly  susceptible  to  matrix  cracking  and  delaminaiions. 

The  modeling  of  laminated  composite  plates  and  shells  has  been  the  subject  of  intensive 
research  in  the  last  two  decades.  Following  the  pioneering  developments  of  Reissner,1,2 
Hildebrand  et  al.,3  and  Mindlin4  in  the  analytic  treatment  of  homogeneous  elastic  plates,  a 
great  number  of  displacement-based,  stress-based,  and  mixed  formulations  for  application  to 
laminate  composites  have  been  explored  (e.g.,  refer  to  the  recent  review  papers  by  Reissner,5 
Reddy,6  Noor  and  Burton,7  and  references  therein).  Reddy6  groups  laminate  theories  into 
three  general  classes:  (1)  equivalent  single-layer  theories  (two-dimensional);  (2)  layer-wise  theo¬ 
ries  (two-dimensional);  and  (3)  continuum-based  theories  (two-  and  three-dimensional).  Of 
the  three  classes,  the  single-layer  displacement  theories  are  the  simplest  and  most  economical 
to  use. 

One  commonly  recognized  drawback  of  the  single-layer  displacement  theories,  however,  is 
that  all  six  components  of  stress,  obtained  from  constitutive  relations,  are  discontinuous  at  the 
ply  interfaces  whereas,  according  to  elasticity  theory,  only  the  inplane  stresses  arc  discontinu¬ 
ous  and  the  transverse  stresses  maintain  continuity  across  the  laminate  thickness.  However, 
when  the  transverse  stresses  are  obtained  by  integrating  the  elasticity  equations  of  equilib¬ 
rium,8  accurate  stress  distributions  can  be  recovered.  To  resolve  the  issue  of  transverse 
stresses  in  a  direct  fashion,  Reissner9,10  proposed  a  mixed  variational  principle  which  uses  the 
three  displacement  components  and  three  transverse  stresses  as  the  independent  variables. 
Although  from  an  analytic  standpoint  this  approach  appears  to  have  some  qualitative  advan¬ 
tages;  from  a  computational  perspective  it  possesses  the  characteristics  deficiency  of  all  stress 
and  mixed  formulations11  which  has  a  relatively  large  number  of  stress  parameters  which  need 
to  be  solved  in  order  to  obtain  element  stiffnesses  with  the  implication  of  an  additional  and 
often  significant  computational  cost.  Thus,  in  large-scale  applications  and,  particularly,  in  com¬ 
putationally  intensive  nonlinear  analyses,  the  elements  of  choice  are  those  that  provide  the 
best  compromise  between  accuracy  and  computational  cost,  with  the  displacement-based 
theories  emerging  as  the  preferred  framework.  In  what  follows  we  narrow  our  focus  on  the 
single-layer  displacement  theories  and  propose  a  new  theory  which  is  specifically  formulated 
with  a  view  on  the  computational  aspects  of  thin  and  thick  composite  laminates. 

In  a  single-layer  displacement-based  theory,  the  basic  assumption  is  that  concerning  the 
through-thickness  approximations  of  the  displacement  components.  The  displacement  compo¬ 
nents  are  expanded  across  the  total  laminate  thickness  with  respect  to  the  thickness  coordi¬ 
nate.  The  expansion  coefficients  (or  the  platc/shell  kinematic  variables)  arc  functions  of  the 
inplane  coordinates  (and  time,  in  dynamics).  Commonly,  the  inplane  displacements  aie 
expanded  with  a  polynomial  of  the  same  degree,  m,  whereas  the  transverse  displacement  expan¬ 
sion  may  be  of  the  different  degree,  n.  Thus,  the  notation  (m,  n}  may  be  conveniently  used 
to  distinguish  between  the  various  single-layer  theories. 
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The  simplest  and  most  extensively  explored  approximation  is  {1,  0}  (i.e.,  a  linear  inplane 
displacement  and  a  constant  transverse  displacement,  totaling  live  kinematic  variables)  or  what 
is  often  referred  to  as  the  first-order  shear  deformable  theories.12'14  These  theories,  which 
are  extensions  of  the  Mindlin  theory,4  enforce  zero  transverse  normal  deformation  by  virtue 
of  n  =  0.  When  formulated  from  the  principle  of  virtual  work,  the  resulting  two-dimensional 
variational  principle  requires  only  C°  continuous  kinematic  fields,  thus  providing  a  convenient 
framework  for  developing  simple  and  computationally  efficient  finite  elements. l3‘ 19 

Other  developments  have  focused  on  higher-order  theories;  those  which  include  transverse 
shear  and  disregard  transverse  normal  deformations  {m>l,  0},6  and  those  that  account  for 
both  transverse  shear  and  transverse  normal  deformations  {m  >  i,  n  >  2},  with  the  latter 
class  of  theories  being  less  prevalent.  "  Although,  in  general,  higher-order  theories  provide 
more  accurate  approximations  of  the  laminate  deformations,  strains  and  stresses,  they  have  not 
been  particularly  suited  toward  finite  element  approximations  due  to  the  presence  of  one  or 
more  limitations  such  as:  (1)  incorporating  a  large  number  of  kinematic  variables  requiring 
C°  or  higher  continuity;  (2)  imposing  natural  edge-boundary  conditions  that  involve  such  non- 
classical  quantities  as  higher-order  stress  resultants;  and  (3)  the  inability  to  model  appropriate 
transverse  stress  boundary  conditions  at  the  top  and  bottom  laminate  surfaces.  Other  com¬ 
mon  deficiencies  include  the  requirement  of  “shear  correction”  factors  that  tune  the  trans- 
verse  shearing  properties,  ’  and  some  theories  lack  a  variational  basis. 

Recently,  Tessler24'26  has  developed  a  higher-order  {1,  2}  theory  for  homogeneous  plates 
which  incorporates  “field-consistent”  transverse  strains  and  is  devoid  of  all  aforementioned  limi¬ 
tations.  The  novel  feature  of  that  theory  is  a  displacement  variational  principle  requiring 
only  C°  and  C'1  continuity  for  the  plate  kinematic  variables,  which  allows  the  development  of 
efficient  plate  finite  elements  having  an  expanded  applicability  range.  One  such  element,  a 
three-node  triangle  with  five  engineering  degrees-of-freedom  (dof)  at  each  node,24  has  demon¬ 
strated  the  same  computational  efficiency  as  its  Mindlin  counterpart.15  In  this  paper,  an 
extension  of  this  theory  to  laminated  composite  plates  and  the  derivation  of  an  efficient  three- 
node  plate  finite  element  are  presented. 

In  Section  2,  the  development  of  the  {1,  2}  laminate  plate  theory  from  three-dimensional 
elasticity  is  presented.  This  is  achieved  by  expand’ng  the  three  Cartesian  displacements  in 
terms  of  the  thickness  coordinate  using  linear  and  parabolic  distributions  for  the  inplane  and 
transverse  displacements.  Additionally,  independent  expansions  are  used  for  the  transverse 
strains.  By  requiring  exact  transverse  stress  boundary  conditions  at  the  top  and  bottom  plate 
faces  and  the  independent  transverse  strains  to  have  equivalent  mean  values  to  those  obtained 
from  the  assumed  displacements  directly,  improved  expressions  for  the  transverse  shear  and 
transverse  normal  strains  are  obtained.  Employing  the  three-dimensional  virtual  work  princi¬ 
ple  yields  a  set  of  seven  partial  differential  equations  of  equilibrium  and  exclusively  Poison 
boundary  conditions.  The  theory  reduces  to  one  of  coupled  lOth-ordcr  stretching-bending  and 
Oth-orrier  transverse  stretching. 

In  Section  3,  a  three-node  plate  element  based  on  the  displacement  variational  statement 
of  Section  2  is  developed.  There  are  seven  kinematic  variables  in  the  formulation;  however, 
because  the  two  higher-order  displacement  variables  do  not  have  spatial  gradients  in  the  varia¬ 
tional  statement,  they  are  -ssumed  to  be  uniform  within  the  element  domain  (C*1  continuous) 
and  are  statically  condensed  out  at  the  element  level.  The  inplane  membrane  displacements 
arc  assumed  to  vary  linearly  across  the  element  while  the  bending  variables  are  interpolated 
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using  an  isoparametric  shape  functions.15  Each  element  node  has  five  engineering  dof 
involving  three  displacements  and  two  normal  rotations.  The  transverse  shear  relaxation 
parameter16  is  also  employed  to  completely  eliminate  shear  locking  in  the  thin  regime. 

In  Section  4,  there  is  a  brief  discussion  on  the  computation  of  accurate  throi’gh-thickness 
distributions  of  the  transverse  stresses  via  a  global  smoothing  technique.  '  In  Section  5. 
analytic  and  finite  element  solutions  to  the  problem  of  cylindrical  bending  of  a  symmetric 
carbon/epoxy  laminate  are  presented.  Results  are  compared  with  Pagano’s  exact  elasticity  solu¬ 
tion.29  In  Section  6,  concluding  remarks  regarding  the  merits  of  the  present  composite  plate 
analysis  are  presented. 


2.  ANALYTIC  BASIS  OF  HIGHER-ORDER  PLATE  THEORY 

The  present  theory  is  developed  in  the  following  manner:  the  Cartesian  displacement 
components  Uj  (i  =  x,  y,  z)  are  expanded  with  respect  to  the  dimensionless  thickness  coordi¬ 
nate  £  =  z/he[-l,l],  where  uz  has  the  special  parabolic  form: 

ux(x,y,z)  =  u(x,y)  +  h59y(x,y)  (a) 

u  (x,y,z)  =  v(x,y)  +  h£6  (x,y)  (b)  (1) 

y  x 

uz(x,y,z)  =  w(x,y)  +  5wx(x,y)  +  (£2+C)w2(x,y)  (c) 

where  £  =  0  designates  the  position  of  the  middle  surface  Sm,  2h  is  the  total  plate  thickness 
and  C  is  a  constant  which  makes  w(x,  y)  a  weighted-average  transverse  displacement  defined 
in  accordance  with  Rcissner,2  that  is, 


h 

“  ■  Sh  |  uz(1'5!>  dz 
-h 


(2) 


Substituting  Equation  lc  into  Equation  2  and  solving  for  C  results  in* 
C  =  -1/5 


(2a) 


The  expansion  coefficients  of  the  inplane  displacements  u,  v,  0y  arc  also  defined  as 
weighted-average  kinematic  variables: 


u 


dz , 


v 


dz 


9  = 
x 


3 

2h3 


u  z 

y 


dz , 


h 

u  z  dz 

J  x 

-h 


(3) 


(A) 


•The  present  nninlion  Inr  w,  dillcrs  from  ituil  in  Reference  24  by  a  factor  of  2/3 
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where  u  and  v  denote  the  midplane  displacements  along  the  x  and  y  directions,  respectively, 
with  9X  and  6y  denoting  rotations  of  a  transverse  normal  about  the  x  and  y  axes,  respectively 
(see  Figure  1). 


H 


X,  u 

X 


Figure  1.  Plate  notation 

For  a  general  composite  layup  composed  of  N  plies,  the  stress-strain  relation  for  each  indi¬ 
vidual  kth  ply  (k  =  1,  2,...,  N)  is  governed  by  a  three-dimensional  Hooke’s  Law  of  the  mono¬ 
clinic  form: 


where  the  elastic  constants  corresponding  to  the  laminate  x-y  coordinates  arc  related  to 

the  ply  principal  direction  constants  as 


C  (k) 
ij 


a .  a .  C 


(k) 


ira  jn  mn 


(6) 
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with  «jj  denoting  the  appropriate  cosine  functions  of  the  coordinate  rotation  between  the  ply 
principal  material  directions  and  laminate  coordinates.30 

The  inplane  strain-displacement  relations  arc  obtained  in  the  usual  manner  as. 


{e,e>Y}  =  {u  ,u  ,  u  -•  u  } 
xx’  yy’  xy  x,x  y,  y  x,y  y,x 


=  {e  +  h 5  k:  ,  c  +  h?  k  ,  v  .  +  he  *  „  } 


y0 


y0  xyQ 


xy0 


(7a) 


=  {u  +  he  8  ,  V  +  he  9  ,  U  +  V  +  he  (9  +0  )} 

x  s  y,  x  ,  y  x  ,y  ,y  ,x  s  x,x  y,y 


A  major  departure  from  a  conventional  displacement  formulation  is  the  way  the  transverse 
strains  are  introduced  into  the  theory.  Here,  the  transverse  strains  arc  expanded  indepen¬ 
dently  in  a  field-consistent  polynomial  form,-4. 


3 


e 

zz 


l 


n=0 


5n  b. 

in 


(i=x,y) 


(7b) 


The  expansion  coefficients  an  and  bj„  arc  obtained  by  requiring  the  weak  transverse  strain 
compatibility 


minimize 


-.2 

u  dz 
z,z 


minimize 


u 

z,  1 


i,z  J 


2 

dz 


(i=x,y) 


(8a) 


(8b) 


*Thc  approach  of  expanding  Ihc  displacement  gradients  u,.,  (i  -  x.  v,  /.)  introduced  in  Reference  24  is  equivalent  to  the  present  one 
based  on  the  transverse  strain  expansions. 
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and  the  satisfaction  of  the  exact  stress  boundary  conditions  at  the  top  (k  =  N)  and  bottom 
(k  =  1)  plate  faces 


x^(x,y,-h)  =  T^(x,y,h)  =  0  (i=x,y) 


(9a) 

(9b) 


where  the  minimization  in  Equation  8  is  performed  witn  respect  to  the  unknown  expansion 
coefficients.  While  Equation  9a  is  the  statement  of  zero  shear  tractions  at  the  top/bottom 
plate  faces.  Equation  9b  results  from  the  transverse  equilibrium  equation  of  three-dimensional 
elasticity 


T(k)  +  T(k) 

xz,x  yz.y 


+  0(k)  =0 
zz  ,z 


(ignoring  the  body  force) 


(10) 


in  which  conditions  Equation  9a  are  incorporated. 

The  resulting  transverse  shear  strains  have  the  parabolic  distribution  satisfying  traction- 
free  boundary  conditions 


■  I  (1-  5!>  KiIo  (I=x,y) 


}  =  {w  +  8  ,  w  +  0  } 
xz0  yz0  ,x  y’  ,  y  xJ 


1  ? 

which  agree  with  Reissner’s  transverse  shear  strains.  ' 

The  resulting  cubic  transverse  normal  strain  can  be  written  as 


(11) 


g  -  e  +  k  *  , 
zz  z0  * 


(12) 


where 


g  =  wx/h 


-  (<  »  <„  ,  K  ,  K  )  =  (Q  ,0  ,  w2/h2,  0  +0  ) 

xo  yo  z0  xy0  y,x  x,y’  2  x,x  y,y' 


^  ( 5  )  ~  Uli  4  2  »  't*  3  > 
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in  which 


<P(Oi 

fil'd 

ii 

(P2U) 

r .  * 
l 

■  gf  [Px(5)/14  +  p}(01  S.} 

r .  = 

l 

(cH/c33)(k: 

=D 

-  (c3i/c33)(k=N) 

s .  = 

i 

(C3./ 

c33)(k: 

=d 

+  (C3./C33)(V'=N)  (i=l ,2,6) 

For  i 

=  3, 

t(Oa 

28 

85 

h  [6Pj 

(O 

-  Pj(5)] 

where  P,  (£)  denote  Lengendre  polynomials: 


PiCO  =  5,  P2CO  -  W-DIl,  P3(0  =  £(5£2-3)/2 


It  can  readily  be  verified  that  for  a  homogeneous  plate  Equations  11  and  12  satisfy  the  three- 
dimensional  eqi"',:hrium  Equation  10  exactly. 

The  plate  equations  of  equilibrium  together  with  the  natural  boundary  conditions  are 
obtained  from  the  three-dimensional  statement  of  virtual  work*  which,  neglecting  body  forces, 
may  be  written  as: 


Tf  /_  (k),  ,  _00.  ,  _(k).  (k).  (kh  (k) 

to  o  e  +  o  6c  +  a  5  e  +  x  <5  y  i-  ",  5  y  +  x  <5y  )  dV 

JJJ  xx  xx  yy  yy  zz  zz  xy  'xy  xz  'xz  yz  'yz 


q*  6u  (x,y,h)  dxdy  + 

Z 


JJ 


q"  <5'i  (x,y,-h)  dxdy 


j  (T  6u  +  T  6u  +  T  5u  )  dsdz  =  0 

!J  x  x  y  y  z  z 


(13) 


Mu'  statement  i  .iii  lx'  retarded  as  a  ’weak"  form  of  the  virtual  work  principle  due  In  the  use  the  mean  conipulil'iluy 
requirement  ter  the  transverse  strains  (see  I  qua! ion  H). 
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where  <5  is  the  variational  operator;  S+  and  S'  denote  the  top  and  bottom  plate  surfaces 
which  are  taken  free  of  shear  tractions  and  loaded  by  normal  pressures, 


T^(x,y,h)  =  0  ( i=x,  y) ,  c/^(x,y,h)  =  q*(x,y)  on  S* 

jlZ  Z  Z 

xiz^(x,y,"h)=  0  (i=x*y)»  °zz^x,y,'h^=  <r(x*y)  on  s" 


(14) 


where  Tj  (i  =  x,  y,  z)  denote  the  prescribed  tractions  on  the  part  of  the  edge  boundary  $ a 
(henceforth,  the  barred  symbols  will  denote  quantities  prescribed  at  the  plate  edges).  Integrating 
Equation  13  across  the  plate  thickness  results  in  the  two-dimensional  plate  virtual  work  principle 


j  j  {  N  5u  +  N  6v  +  N  (6u  +  5v  )  +  N  5(w-,/h) 

JJ  [x,x  y.y  xy  ,  y  ,  x  z  1 


in 


+  M  60  +  M  69  +  M  (60  +60  )  +  M  6(w7/h2) 

x  y,x  y  x,y  xy  x,x  y.y  z 


+  Q  (6w  +  60  )  +  Q  (6w  +  50  )  1  dxdy 

x  ,  x  y  y  •  y  x  J 

-  JJ  jqi(<5w  +  |  6w2)  +  q26w1  |  dxdy 


m 


-f  {  Sxn 


6u  +  N  6v  +  M  60  +  M  50  +  Q  6w  +  Q  6w,  +  Q  5w, >  ds  =  0 

yn  m  \m  v  rn  t  a  _  *  I 


yn 


xn  y  yn  x  zn 


(15) 


where  the  plate  stress  resultants  are  defined  in  Appendix  A.  Integrating  Equation  15  by 
parts  results  in 

If  f  -[N  +  N  ]  6u  -  [N  +  N  ]  6v 

JJ  [  x.x  xy,y  y.y  xy ,  x 

S 

ra 

+  [  -  M  -M  +  Q  ]  50  +  [-  M  -M  +  Q  ]  60 
x , x  xy,y  ^x  y  y,y  xy,x  y  x 


[Qx>x+  Qy,y+  ^  5W  +  fN7/h  '  ^  5wl 


+  [iWh2  -  g  q x ]  6w2  j  dxdy 
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+  <j>  I  (N  -  N  )  5u  +  (N  -  N  )  6v  +  (M  -  M  )  59  +  (M  -  M  )  59 
J  [  xn  xn  yn  yn  xn  xn  y  yn  yn  x 

C 

o 

+  (Q  -  Q  )  6w  +  Q  6w,  +  Q  6w2l  ds 

zn  zn  zx  z2  J 


HN  5u  +  N  5v  +  M  59  +  M  59  +  Q  5w  1  ds  =0 

xn  yn  xn  y  yn  x  zn  J 


(16) 


where 


N  = 

N  n 

+ 

N  n  , 

N 

=  N  n 

+ 

N  n  , 

xn 

X  X 

xy  y 

yn 

xy  x 

y  y 

M  = 

M  n 

+ 

M  n  , 

M 

=  M  n 

+ 

M  n 

xn 

X  X 

xy  y 

yn 

xy  x 

y  y 

and 


(17) 


(n  ,  n  }  =  (cos(x,n) ,  cos(y,n)} 
x  y 


where  the  Ca  and  Cu  are  parts  of  the  boundary  C  (QUC^  =  C  and  CunCff  =  0  )  surround¬ 
ing  the  middle  surface  Sm  where  the  tractions  and  displacements  arc  prescribed;  « nd  n  and  s 
denoting  the  outward  normal  and  tangential  coordinates. 

Thus,  the  variational  principle  provides  seven  equations  of  equilibrium  written  as: 


(6u) : 

N  +  N 
x,x  xy, y 

*  0 

(18a) 

(6v): 

N  +  N 

xy,x  y.  y 

=  0 

(18b) 

(S9y): 

M  +  M  -  Q 

x,x  xy , y  x 

=  0 

(18c) 

(6V: 

M  +  M  -  Q 

xy,x  y, y  My 

=  0 

(18d) 

(6w) : 

Qx,x  +  Qy,y  +  qi 

=  0 

(18e) 

(6wx) 

Nz/h  -  q2 

=  0 

( 1 8  f ) 

(<5w2): 

Mz/h2  -  l  qx 

=  0 

(18g) 
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and  the  associate  boundary  conditions: 


fa')  Poisson  boundary  conditions 


N  =  N  ,  N  =  N  ,  M  =  M  ,  M  =  M  ,  0=0  onC 

xn  xn  yn  yn  xn  xn’  yn  yn*  wzn  wzn  on  a 


(19) 


(h)  Vanishing  higher-order  shear  resultants 


q  =  q  s  o  on  C 

2  ^  2  2  cr 


(20) 


(C-)-Displacement  boundary  conditions 


u  =  u,  v  =  v, 


on  C 

u 


(21) 


The  plate  constitutive  equations  relating  stress  resultants  to  strains  are  summarized  in  Appen¬ 
dix  B.  Upon  their  introduction  into  Equation  18  there  result  two  sets  of  equilibrium  equa¬ 
tions  in  terms  of  displacement  variables: 


qA2u  ■  Lalq  (22) 

"  =  RbLblu  +  V  (23) 


where 


T 

q 


q2). 


T 

u 


{u,  V,  w. 


9  }  and 

y 


Wi.  w2> 


(24) 


The  Qa,  Rh.  and  Sn  are  various  stiffness  coefficient  matrices.  The  Laj  and  Lm  are  first-order 
linear  differential  operators,  and  the  La2  is  a  second-order  linear  differential  operator.  The 
explicit  forms  of  these  matrix  operators  arc  rather  cumbersome  and,  for  this  reason,  they  are 
not  presented  here.  Note  that  Equation  22  contains  five  second-order  coupled  partial  differen¬ 
tial  equations  of  equilibrium  in  the  five  displacement  functions  u,  v,  w,  dx  and  6y  (i.c.,  a  l()th- 
order  system),  whereas  Equation  23  has  two  zero-order  equations  in  the  w |  and  w2 
displacements. 
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The  solution  to  a  typical  plate  boundary-value  problem  is  obtained  by  solving  Equation  22 
for  u  subject  to  the  boundary  conditions  of  Equations  19  and  21,  and  then  solving  the  “auxil¬ 
iary”  system  of  Equation  23  for  w.  The  computation  of  stress  resultants  follows  from  the 
plate  constitutive  equations  (Appendix  B),  while  the  stresses  are  computed  from  the  three- 
dimensional  Hooke’s  Law  in  Equation  5.  In  Section  4,  there  is  a  brief  discussion  of  an  alter¬ 
native  procedure  for  recovering  accurate  transverse  shear  stresses. 

2.1  Remark 

The  plate  theory  can  be  reduced  to  a  shear-deformable  one  by  enforcing  infinite  rigidity 
in  the  transverse  normal  material  direction  (C33  =  «)  and  setting  =  C23  =  0.  Further¬ 
more,  with  the  additional  assumption  of  infinite  rigidity  in  the  transverse  shear  directions  (C44 
=  C55  =  oo),  the  theory  yields  results  of  the  classical  laminate  plate  theory  (see  Section  5). 

3.  FINITE  ELEMENT  DEVELOPMENT 

The  particular  form  of  the  variational  principle  (see  Equation  15)  lends  itself  well  to  the 
development  of  simple  and  efficient  finite  elements  comparable  in  computational  efficiency  to 
the  first-order  shear-deformable  elements.15*18  Recognizing  that  Equation  15  contains  gradi¬ 
ents  of  u,  v,  w,  #x  and  6y  that  do  not  exceed  order  one  implies  C°  continuity  of  these  vari¬ 
ables.  In  addition.  Equation  15  has  no  gradients  of  wj  and  w2;  therefore,  C’1  continuous 
approximations,  which  arc  discontinuous  at  the  element  boundaries,  can  be  employed  for  these 
variables. 

An  important  simplification  due  to  wj  and  W2  being  C"1  continuous  is  that  their  dof  can 
be  conveniently  eliminated  at  the  element  level  via  static  condensation.  Alternatively,  the 
elimination  of  wj  and  w2  can  be  performed  in  the  variational  statement  by  realizing  that  in 
Equation  15  the  grouped  terms  associated  with  the  arbitrary  variations  c5wj  and  dw2  must  van¬ 
ish  independently.  The  result  is  Equation  18f  and  Equation  18g  which,  by  using  the 
constitutive  relations  for  Nz  and  Mz  (Appendix  B),  can  be  readily  solved  for  w!  and  w2;  these 
are  then  replaced  into  Equation  15  prior  to  the  finite  element  approximation.  With  this 
approach,  Equation  18f  and  Equation  18g  are  satisfied  exactly  whereas  the  static  condensation 
at  the  finite  element  level  accomplishes  the  same  purpose  consistent  with  the  finite  element 
approximations. 

3.1  Thin-Regime  Considerations 

By  expressing  Equation  15  in  dimensionless  form,  it  can  further  be  shown  that  this  varia¬ 
tional  principle  is  of  a  penalty-constraint  type  analogous  to  that  of  the  Mindlin  theory  (e.  g., 
refer  to  Reference  18).  In  the  thin  regime  as  h  -»  0,  it  enforces  the  Kirchhoff  constraint  of 
the  vanishing  transverse  shear  angles,  that  is 


(w  +  0>  w  +  8  }  -►  0 
,  x  y*  ,y  x' 


(25) 


and  the  inextensibility  of  the  transverse  normal 


11 


(26) 


Wlf  w2  -*■  0 

Whereas  Equation  26  imposes  no  restrictions  upon  element  interpolations  for  wj  and  W2,  the 
Kirchhoff  constraints  (see  Equation  25)  require  field  consistency  in  the  interpolations  for  w, 
dx  and  6y  which  may  be  achieved  through  the  use  of  anisoparametric  shape  functions.15'18  In 
addition,  an  element-level  relaxation  of  the  transverse  shear  penalty  parameter  may  be  neces¬ 
sary  when  the  field  consistency  is  violated;  for  example,  when  excessive  kinematic  boundary 
restraints  exist.16 

In  what  follows,  anisoparametric  interpolations  for  a  three-node  element  are  summarized, 
and  a  brief  discussion  on  the  implementation  of  shear  relaxation  is  presented.  To  distinguish 
element  approximations  from  their  analytic  counterparts,  the  former  are  superscribed  with  an  l 
which,  in  the  present  notation,  signifies  a  characteristic  length  scale  of  the  discretization. 

3.1.1  Three-Node  Anisoparametric  Interpolations 


i,  ,  ;  r  i 

u  (x,y)  =  E  t.u. , 

i=l 


e*  <x.y) 


3 

=  I 
i=l 


.0  . 
i  xi 


v  (x,y) 


3 

=  I 
i=l 


,  1 
C  .V. 
1  1 


ej  <x,y) 


3 

=  s 

i=l 


Cfl  Linear 

Shape 

Functions 


(27) 


where  £i  are  the  area-parametric  coordinates  given  as 

L.  =  (c,  +  b.x  +  a.y)/2A  (A  is  the  triangle  area) 
11  1  1 


and  ai,  bi,  and  Ci,  are  coefficients  of  nodal  coordinates 
ai  =  W  ci  =  xjWj 


^(x.y)  -  +  e;.(Vrt,  -  bjNk+3)/8 


+  9yi(aj\+3  -  VW/81 


C°  Quadratic 

Shape 

Function 


(28) 


where  Nj  are  quadratic  shape  functions 


N.  =  C.(2C.  -  1),  N  _  =  4C.C. 
i  l  i  i+3  i  j 
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and  the  subscripts  assume  the  values: 


i=  1,2,3;  j=  2,3,1;  k-3,1,2. 


(x,y)  =  W* 
(x,y)  =  Wj 


C'1  Uniform 
Functions 


(29) 


Introducing  Equations  27  through  29  into  Equation  15,  and  performing  the  necessary  varia¬ 
tions  with  respect  to  the  displacement  dof  and  integrating  over  the  element  domain,  the  ele¬ 
ment  stiffness  equilibrium  equations  are  derived.  The  and  dof  are  then  eliminated 
using  static  condensation  yielding  a  three-node  clement  with  five  engineering  dof  at  each 
node;  i.e.,  {uj,  Vj,  Wj,  0^,  0yi}  (i  =  1,  2,  3)  (see  Figure  2). 


KEY:  o  u,  i,  0^,  9^  do!  □  t  |  dof 


Figure  2.  A  three-node  stretching-bending  plate  element  with 
transverse  shear  and  transverse  normal  deformation. 


3.1.2  Shear  Relaxation 

The  shear  relaxation  parameter,  (p  s,  is  a  multiplier  introduced  in  the  transverse  shear  con¬ 
stitutive  equations  for  the  element 


Q£  -  G  y2. 


(30) 


The  parameter  is  determined  from  the  diagonal  transverse  shear.  ks0,  and  bending,  k^,  stiff¬ 
ness  coefficients  corresponding  to  rotational  dof  of  the  "unrelaxed"  element  stillness  matrix:15'18 
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1 


(31) 


4> 


2 

S 


1  +  2  (S  k9/Z  k9) 
s  b 


where  the  summations  are  carried  over  all  six  diagonal  terms.  The  use  of  shear  relaxation 
avoids  locking  in  the  thin  regime. 


The  computer  implementation  of  this  three-node  plate  element,  enabling  its  use  with  the 
general  purpose  finite  element  code  ABAQUS,  is  discussed  in  Reference  31. 


4.  TRANSVERSE  STRESS  CALCULATIONS 

Due  to  their  heterogeneous  nature,  laminated  composite  plates  generally  exhibit  ply- 
interface  discontinuities  in  the  inplane  stresses,  as  well  as  in  the  transverse  strains.  Also,  the 
inplane  strains  ire  continuous  across  the  laminate  thickness  and  the  transverse  stresses  arc 
also  continuous,  although  their  through-thickness  gradients  may  be  discontinuous  at  the  ply 
interfaces. 

Since  the  present  theory  relies  upon  continuous  strain  distributions  across  the  entire  lami¬ 
nate  (as  do  all  single-layer  displacement  theories),  accurate  inplane  stresses  can  be  recovered 
from  the  ply  constitutive  law  (see  Equation  5);  however,  the  direct  recovery  of  transverse 
stresses  in  this  manner  are  often  erroneous  producing  nonphysical  discontinuities  at  the  ply 
interfaces. 

To  obtain  accurate,  through-thickness  continuous  transverse  shear  stresses,  one  needs  to 
integrate  the  three-dimensional  equilibrium  equations  using  the  inplane  stresses  obtained  from 
the  constitutive  relations  (see  Equation  5  and  refer  to  Reference  8),  that  is 


„(k>> .  r^)  +  T<">  >  dz 

xz  J  xx, x  xy,y 


^12»  C13,  C16^  — 


XX 


yy 


zz 


xy 


+  (C16,  C j g ,  C36,  C66^  ^  a 


3y 


XX 


yy 


zz 


xyj 


dz. 


(31a) 
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where  the  strains  are  defined  in  Equation  7a  and  Equation  12.  Having  obtained  the  trans¬ 
verse  shear  stresses,  the  transverse  normal  stress,  o^\  can  be  found  either  by  integrating 
Equation  10  or  applying  Hooke’s  Law  (sec  Equation  5).  As  will  be  shown  in  Section  5,  in 
this  plate  theory  both  approaches  yield  accurate  predictions  for  o72^\  however,  the  Hooke’s 
Law  stress  recovery  yields  slightly  discontinuous  transverse  normal  stresses  at  the  ply  interfaces. 

o 

While  in  the  context  of  this  plate  theory,  as  well  as  others,  the  equilibrium  stress-recov¬ 
ery  approach  has  been  effective;  its  application  to  the  finite  element  method  is  not  straightfor¬ 
ward.  For  instance,  in  the  case  of  the  present  finite  element,  all  strain  gradients  in  Equation 
31  vanish  identically  (the  strains  are  constant  within  the  element  domain);  therefore,  the 
direct  utilization  of  these  finite  element  strain  gradients  will  always  result  in  zero  transverse 
shear  stresses. 

To  overcome  this  difficulty,  a  global  smoothing  procedure  ’  yielding  continuous  solu¬ 
tions  for  the  strains  and  their  normal  gradients  across  the  entire  plate  domain  was  employed. 
The  strain  gradients  were  then  used  in  Equation  31  to  compute  accurate  through-thickness  dis¬ 
tributions  of  the  transverse  shear  stresses.  For  complete  details  on  the  global  smoothing  pro¬ 
cedure  the  reader  is  referred  to  Reference  28. 

5.  RESULTS  AND  DISCUSSION 

Numerical  studies  are  carried  out  for  the  problem  of  cylindrical  bending  of  a 
carbon/epoxy  symmetric  angle-ply  ([30/-30]s)  laminate  subjected  to  a  sinusoidal  transverse  pres¬ 
sure  q  =  qosin(jrx/L)  for  which  an  exact  elasticity  solution  is  available29  (see  Figure  3a). 

This  problem  is  rather  challenging,  exhibiting  significant  inplane  shear  coupling  in  each  ply  of 
the  laminate.  The  ply  material  properties  arc  taken  as 
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(32) 


EL  =  25  x  106  psi,  Et  =  106  psi 

Glt  =  0-5  x  106  psi,  Gjt  =  0.2  x  106  psi 


i>LT  —  i^xT  =  0.25 


where  L  and  T  denote  the  longitudinal  and  transverse  ply  material  directions,  respectively. 

The  results  from  the  study  were  obtained  analytically  from  the  present  Higher-Order 
Uieory,  shown  as  HOT-ANALYTIC  (for  the  analytic  solution  approach  refer  to  Reference 
24),  an  exact  elasticity  approach,29  designated  as  EXACT,  the  classical  laminate  plate  theory 
(CPT),  and  the  present  finite  element  analysis,  designated  as  HOT-FEM.  Also,  results  corre¬ 
sponding  to  the  shear-deformable  version  of  the  present  theory  (see  Remarks  2.1)  are  pre¬ 
sented  and  labeled  in  the  figures  as  SHEAR-DEFORMABLE. 

Figure  3b  depicts  a  twenty-element  finite  element  discretization  of  a  narrow  strip  spanning 
one  quarter  of  the  loading  wavelength  (i.e.,  L/2).  Appropriate  symmetric  boundary  conditions 
and  multiple-point  kinematic  constraint  equations  are  enforced  to  simulate  the  cylindrical  bend¬ 
ing  of  the  plate.  In  Figure  4,  the  percent  error  of  the  maximum  midplane  transverse  dis¬ 
placement  is  plotted  versus  the  half  wavelength  of  loading-to-thickness  ratio  (L/2h).  As  L/2h 
decreases,  the  observed  deflection  becomes  dominated  by  transverse  shear  and  transverse  nor¬ 
mal  deformations.  The  thick  regime  improvement  in  accuracy  afforded  by  the  inclusion  of 
both  of  these  effects  in  the  present  tneory  is  clearly  evident. 

The  remainder  of  the  discussion  concerns  the  strain  and  stress  predictions  for  the 
moderately  thick  (L/2h  =  10)  and  thick  (L/2h  =  4)  deformation  regimes.  Figures  5  and  6 
depict  through-thickness  distributions  of  the  maximum  inplane  (fH  and  y^)  and  transverse  yxz 
and  ezz)  strains,  respectively.  The  observed  departures  from  the  exact  elasticity  solutions, 
particularly  at  the  outer  fibers  for  L/2h  =  4,  can  be  attributed  to  the  linear  distributions  for 
ux  and  uy  in  the  present  theory.  Note  that  their  thickness  variations  according  to  the  elastic¬ 
ity  solution  become  increasingly  nonlinear  in  the  thick  regime. 

Figure  7  shows  the  distributions  of  the  maximum  transverse  shear  txz  and  transverse  nor¬ 
mal  Ojz  stresses.  These  stresses  were  computed  following  the  conventional  Hooke’s  Law 
stress  recovery  (see  Equation  5),  designated  as  HOT-HOOKE,  and  by  integrating  the 
three-dimensional  equilibrium  equations  of  elasticity  (sec  Equation  31),  designated  as 
HOT-EQUIL.  It  is  of  interest  to  note  that  while  the  equilibrium-based  transverse  shear 
stress  (rxz)  recovery  yielded  superior  results  over  those  obtained  from  Hooke’s  Law,  the  trans¬ 
verse  normal  stress  (azz)  recovery  is  highly  accurate  by  both  methods. 
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Figure  4.  Percent  error  of  maximum  midplane  deflection  versus  l_/2h  ratio. 
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L/2h  =  10 


L/2h  =  10 


U/q0]Ozz(L/2,y,z) 


.  5  0  *■*=— - - - 

0  •  .20  .40 


[2h/q0L]Txz (L, y, z) 


L/2h  =  4 


Figure  7.  Distributions  of  transverse  stresses  across  thickness  for  L/2h  =  10  and  4 


6.  CONCLUDING  REMARKS 

A  computationally  viable  {1,  2}  higher-order  plate  theory  for  stress  analysis  of  composite 
laminates  was  presented.  The  theory  incorporates  both  transverse  shear  and  transverse  nor¬ 
mal  deformations  and  is  based  on  linear  through-thickness  expansions  for  the  inplane  displace¬ 
ments  and  a  special  parabolic  distribution  for  the  transverse  displacement.  In  addition, 
independent  expansions  were  used  for  the  transverse  strains,  satisfying  exact  transverse  stress 
equilibrium  at  the  top  and  bottom  plate  surfaces  and  having  equivalent  mean  values  to  those 
obtained  directly  from  the  assumed  displacements.  From  the  three-dimensional  virtual  work 
principle,  the  two-dimensional  plate  variational  principle  is  derived,  giving  rise  to  a  coupled 
l()th-order  stretching-bending  theory  subject  to  exclusively  Poisson  boundary  conditions  and 
two  “auxiliary”  Oth-ordcr  transverse  stretching  equations. 

To  demonstrate  the  theory’s  computational  aspects,  an  efficient  three-node  stretching- 
bending  plate  element  was  developed.  It  was  shown  that  the  higher-order  transverse  displace¬ 
ments  can  be  eliminated  cither  directly  from  the  variational  statement  or  via  static 
condensation  at  the  element  level.  With  the  use  of  the  latter  approach,  a  5-dof-per-nodc  ele¬ 
ment  was  produced,  having  the  same  computational  efficiency  as  a  comparable  first-order  shear 
deformable  clement.  In  addition,  an  effective  approach  for  the  recovery  of  reliable  transverse 
stresses  was  discussed. 
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In  closing,  the  practical  benefits  of  the  proposed  laminate  plate  analysis  are:  (1)  there 
are  no  requirements  for  “shear  correction”  factors;  (2)  reliable  three-dimensional  distributions 
of  displacements,  strains  and  stresses  are  attainable;  (3)  the  range  of  applicability  includes 
both  thin-  and  thick-section  composite  laminates;  and  (4)  the  finite  element  formulation  repre¬ 
sents  a  viable  enhancement  over  the  widely  used  first-order  shear-deformable  models. 
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APPENDIX  A.  PLATE  STRESS  RESULTANTS  AND  STRAINS 


Stress  Resultants 


N  h. 
r  c  * 


NT={N,N,N,N  }=y  f  {a  ,o  ,  a  ,  a  } 

x’  y*  z’  xy  £  J  xx’  yy’  zz’  xy 


(k) 


dz 


k=l  h, 


k-1 


MT  =  {M  ,  M  ,  M  ,  M  }  =  V  f  {(za(k)+  *1a(k)),  (zo(k)+  $2a(k)), 
x’  y*  z*  xy  l  J  xx  zz  yy  y2  zz  ’ 

k=1  Vl 


_(k)x  /  (k)  (kh. 

(<j>3azz  ),  (zxxy  +  4>6azz  )}  dz 


QT  -  (Qx.  Qy)  -[  {h|‘  |  (1 

k“  Vl 


£2){t(1c),  T(k)}  dz 
xz  yz 


Plate  Strains 


e  =  {e  ,  e  ,  e  ,  y  }  =  {u  ,  v  ,  w,/h,  u  +  v  } 

x0  y0  z0  xy0  ,  x*  ,  y  1  ,  y  ,xJ 

T 

k  =  {k  ,  k  ,  k;  ,  tc  }  =  {e  ,  0  ,  w2/h2,  9  +0  } 

X0  yo  z0  xy0  y.x  x,y  2  x,x  y,y' 


f  -  (y  >  Y  }  =  {w+0,  w+8} 

1  xz  q ’  yz0  , x  y’  ,y  x' 


(A.  1 ) 


(A. 2) 


(A. 3) 

(A.  4) 
(A. 5) 
(A. 6) 


^oplidU  Normal  Irasliims 

q i  =  q*-  q'»  q2  ■  q*+  q 


(A. 7) 
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(A.o; 


r 


h 

{  <v  y 2  dz 

-h 


(A. 9) 


h 

■ 

{Q  .  6  ,  6  }  =  T  {1,  c,  C2  -  1/5}  dz 

zn  z  i  z  j  J  z 

-h 


(A.  10) 


Note:  From  (A.  10),  one  can  readily  verify  that  the  vanishing  conditions  on  QZ1 
(i  =  1,  2)  in  Equation  20  imply  that  Tz  must  vary  paraholically  as  (1  -  t2)  across  the  plate 
thickness. 


APPENDIX  B.  PLATE  CONSTITUTIVE  RELATIONS 


Invoking  Hooke’s  Law  (see  Equation  5)  into  the  plate  stress  resultants  (Appendix  A) 
yields  the  plate  constitutive  relations  which  can  be  expressed  in  matrix  form  as 


r  hi 

A  BO 

E 

M 

= 

T 

B  D  0 

1C 

.  Q  . 

0  0  G 

.  T  . 

(B.l) 


where  the  vectors  of  plate  stress  resultants  and  strains  are  defined  in  Appendix  A;  the 
constitutive  matrices  of  elastic  coefficients  are: 

Membrane  Elastic  Stiffness 


A  -  [A  ]  (i,j  =1,2, 3, 6) 


where 

K 

f  C(k>  dZ 

J  1J 


N 


k=l 


(3.2) 


Membrane-Bending  Coupling  Elastic  Stiffness 


B  =  [B_]  (i,j  =  1,2, 3, 6),  (B.3) 

where 

N  \ 

Biq  =  L  j  {  Clq>z  +  Ci5\  ]  dz  =  1.2, 3, 6!  q-  1,2,6) 

k"1  \-l 

N 

Bi3  =  J  f  Ci3^3  dz  Ci  =  1.2, 3, 6) 

k=1  \-l 
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Bending  Elastic  Stiffness 


D  =  [D..]  (i,j  =  1,2, 3, 6), 
where 

N 

d  -y 

qq  L 

k=l 


“33-1 
k=1  Vi 


v^N  he- 

k=1  Vi 


N  ^ 

V  =  I  I  {  Cq3)*3Z  +  C33}V3}  dZ  (q  =  2*  6) 

k=1  Vl 


(C(k)<j>  +  C(k)c 

q3  Yr  r3 


i  )z  +  C^k)< 

q  33 


q  r 


(q,r  = 


n 


V 


c(kV 

qq 


+  2  C 


00, 

q3 


i  z  +  c$kV 
q  33  Yq 


dz  for 


Transverse  Shear  Elastic  Stiffness 
G  =  [G.j]  (i,j  =  5, A), 

where 

Gij'I  }l'Ci5)[-T<1  - 

k=1  \-l 


(B.A) 


-  1,2,6 


|  dz 

1,2,6) 


(B.5) 
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